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ABSTRACT 


For univariate random sequences, the power spectral density acts like a proba- 
bility density function of the frequencies present in the sequence. This dissertation 
extends that concept to bivariate random sequences. For this purpose, a function 
called the joint spectral density is defined that represents a joint probability weight- 
ing of the frequency content of pairs of random sequences. Given a pair of random 
sequences, the joint spectral density is not uniquely determined in the absence of any 
constraints. Two approaches to constraining the sequences axe suggested: (1) assume 
the sequences axe the margins of some stationary random field, (2) assume the se- 
quences conform to a particular model that is linked to the joint spectral density. For 
both approaches, the properties of the resulting sequences axe investigated in some 
detail, and simulation is used to corroborate theoretical results. It is concluded that 
under either of these two constraints, the joint spectral density can be computed from 
the non-stationary cross-correlation. 
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CHAPTER 1 
INTRODUCTION 


1-1 Problem Statement 

The power spectral density (PSD) of a random sequence can be loosely inter- 
preted as a probability density function (PDF) of the frequency content of that se- 
quence. Stated more precisely, if the PSD of a random sequence *(() is denoted by 
/(A), then /(A)dA represents the expected value of the portion of the total power in 
x(() due to components with frequency in the interval (A, A + dXj [1], Normalizing 
/(A) by the total power gives a probability weighting for each frequency analogous 
to a PDF of the frequency. As such, the power spectral density has all the properties 
of a probability density function. This is one of many parallels between the spectral 
theory of random sequences and the probability theory of random variables. The inte- 
gral of the PSD, sometimes called the integrated spectrum or the'spectral distribution 
function, is analogous to a probability (cumulative) distribution function. Also, in 
the time domain, the auto-covariance is closely related to the class of characterist.c 
functions. Much of the terminology echoes the similarity between concepts in these 
two bodies of theory. In this context it makes sense to talk about the statistics of the 
spectrum where frequency is treated as a random variable. For example, computing 
the spectral moments is often of considerable interest in applications such as radar 
signal processing (2). Covariance based approaches to estimating the spectral mean 
and width can be understood in terms of the moment theorem from probability the- 
ory applied to covariances and spectral densities instead of characteristic functions 
and probability densities [3]. 

This one-to-one correspondence begins to break down when the joint statistics 
of more than one random sequence are under consideration. Consider the bivariate 
case of two random sequences. What we would like to have is a single function of 


two frequency arguments that acts like a joint probability density function ,„ r the 
frequency content in the two sequences. Conventional spectral analysis of bivariate 
random sequences involves four functions of a single frequency argument that are 
typically written as a twoby-two matrix valued function called the spectral density 
matnx. The PSDs of the two sequences, sometimes called their auto-spectra, are 
placed along the majn diagonal and their cross-spectra are placed along the cross di- 
agonal. In essence, the auto-spectra correspond to the marginal probability densities 
of the frequency content for each sequence, and the cross-spectra contribute infor- 
mation about the relationships between the two sequences. It is well known from 
probability theory that the joint probability density function cannot be uniquely de- 
termined based on knowledge of the marginals [4J. Thus the ability to form this 
desired function via conventional bivariate analysis depends on whether the cross- 
spectra contain the extra information that is necessary to uniquely determine the 
joint probability density function. 

This function that acts like a joint probability density function of the frequency 
content in the two sequences will be called the joint spectral density (JSD). The 
terminology, joint spectral density, requires some clarification because the individual 
words (joint, spectral, and density) are rather generic and have been used in various 
combinations to denote things that are not intended here. To help reduce the risk 
of confusion, first consider what is NOT meant by the joint spectral density as the 
term is defined in this dissertation. It is not the spectral density matrix of a bivariate 
random function. It is not the coherency, or the magnitude squared coherency. It is 
not the bispectrum or any other higher-order spectrum. What is meant by the joint 
spectral density is a joint probability weighting of the frequency content in a set of 

two or more random sequences that is completely analogous to a joint probability 
density function. 


Why bother about this joint spectral density function when there exists a well 
developed theory for multivariate random functions based on the covariance matrix 
in the time domain and the spectral density matrix in the frequency domain? The 
short answer to this question is that having the joint spectral density will make all 
the techniques for manipulating PDFs available to the spectral analysis of bivariate 
random sequences. In particular, random variable transformations, or the so-called 
algebra of random variables, are useful because they can be used to determine PDFs 
of functions of the frequency in two random sequences. If x and y are any two random 
variables, which may represent the frequency content in two random sequences, and 
</(•, •) is a function, then the random variable 

h = 

can be expressed in terms of g(-, •) and the joint probability density function of x and 

V [ 4 ]- 

The focus of this dissertation is to form a precise definition of the joint spectral 
density function in terms of the spectral representation of random sequences, explore 
its properties, and determine under what conditions (if any) it can be estimated from 
realizations of a bivariate random sequence. 

1.2 The Radar Problem 

An application where the joint spectral density could be used to great advan- 
tage is in the processing of radar backscatter return from distributed targets such as 
weather. A pulsed Doppler weather radar can be used to construct a map of the radial 
windspeed 1 over a scanned air volume (see Figure 1.1) [2], Each of the boxed regions 
in Figure 1.1 represent a range cell volume from which samples of the backscatter re- 
turn are collected over time (see Figure l .2). These samples form a random sequence 

‘The component of the vector wind on a radial from the radar platform will be referred to as 
the radial windspeed. 
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which the frequency content is related via the Doppler shift to the velocities of 
wind-blown particles within the range cell volume. 

An important part of the signal processing task involves estimating the zeroth, 
first, and second moments of the Doppler spectrum in each range cell [5], These three 
moments of the Doppler spectrum correspond to the signal power, mean velocity, and 
spectrum width respectively. When the goal of signal processing is the detection of 
certain types of weather events (e.g., gust fronts, tornadoes, windshears), the radial 
windspeed gradient, g, is often of particular interest. Knowledge of the joint spectral 
density of two adjacent range cells would make it possible to compute a probability 

density of windspeed gradient at each range increment through the use of random 
variable transformations. 

An example of a function of the windspeed gradient that is of interest in the 
wmdshear detection problem, is a hazard index called the “F”-factor given by 

f = — — 1 4- 

9 dr [ v 3 v a \ ’ 

where V, is the aircraft groundspeed, V. is the aircraft airspeed, h is the altitude. 

and g is the acceleration due to gravity [6J. The radial windspeed gradient can be 
estimated by 

dv 1 

^ = ^( 0 ,- 0 ,), 

where t>, and v, are the windspeeds in two range cells and Ar is the distance between 
them. If the joint probability density of u, and v 2 is known, then it is a simple ma.ie, 
to derive a probability density of hazard factor as the density of the difference of two 
random variables scaled by a constant [4). The JSD provides a way to get that jomi 
probability density of the windspeeds in two range cells. 

1.3 A Note on Terminology 

Much of the terminology associated with the spectral analysis of random .. 
quences is used inconsistently in the various texts and papers available on the - . . 

or---- 

C .i U ■ v fl*. AiU 1 if 



ject. In this section, ground rules are established for how various terms are used in 
this dissertation. First, there is no consistent way of distinguishing between random 
functions of a continuous parameter (sometimes time, but not restricted to be so) 
and random functions of a discrete parameter. In this regard, a convention follow- 
ing that of Yaglom is used so that random functions of a continuous parameter are 
called random processes and random functions of a discrete parameter are called ran- 
dom sequences [7], When the argument may be either continuous or discrete it is 
called a random function. Throughout most of what follows, random sequences are 
treated exclusively, however, in most cases the results are equally applicable to ran- 
dom processes with only minor adjustments (usually involving a change of the limits 
of integration from (-t, t r] to (-00, 00]). 

When more that one random function is observed, they are collectively referred 
to as a multivariate random function. Single random functions of more than one ar- 
gument are sometimes called multidimensional random functions, but here the term 
random field is preferred. The term random field makes no distinction between con- 
tinuous and discrete arguments, but unless otherwise specified, the arguments are 
assumed to be discrete. 

The exact meaning of margin as it is used in this dissertation depends on the 
type of function to which it is applied. For random fields and their correlation func- 
tions, the margin refers to the function that results from setting one or more of the 
arguments equal to zero. The margins of spectral distribution functions result from 
sett.ng one or more of the arguments equal to x (or co in the continuous case). Fi 

nally, for spectral density functions the marg.n is formed by integrating out one „r 
more of the arguments. 


1.4 Contribution to the Field 

The contributions of this dissertation are twofold. First, a definition of the joint 
spectral density is given in terms of the spectral representation of random field*. 
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setting the problem on firm mathematical ground. According to this definition, the 
joint spectral density is not uniquely determined by the two sequences. Second, two 
classes of sequences are identified that enable the joint spectral density to be uniquely 
determined. The class of marginal sequences is considered first. Several theorems are 
offered on the properties of marginal sequences, and examples are provided to illus- 
trate those properties. The second class of sequences is based on an exponential model 
that explicitly incorporates the joint spectral density in its definition. Simulations 
are used in both cases to corroborate the theory. 

1.5 Organization 

In the second chapter, the spectral and correlation theory of univariate random 
sequences is reviewed to establish basic concepts that axe extended to bivariate ran- 
dom sequences and random fields in Chapter 3. Also in Chapter 3 a definition for 
the joint power spectral density is proposed based on establishing a connection be- 
tween the spectral representation of a bivariate random sequence and the spectral 
representation of a random field. In Chapter 4 this connection is made by assuming a 
particular relationship between the bivariate random sequence and a stationary ran- 
dom field. Another approach to constraining the problem is explored in Chapter i 
where the bivariate random sequence is assumed to conform to a special model that 
allows the joint power spectral density to be computed under some special circum 

stances. Finally, in Chapter 6 some conclusions are drawn from this research and 
suggestions axe made for future work. 




CHAPTER 2 

THEORY OF STATIONARY RANDOM SEQUENCES 

In this chapter basic results in the correlation theory of stationary random se- 
quences are reviewed. Sections 2.1 thru 2.3 provide some necessary background mate- 
r,al from probability theory concerning functions that characterize random variables 
and, by extension, sequences of random variables. Section 2.4 is a mathematical 
aside defining the Stieltjes integral which occurs frequently in this and the following 
Chapters. The remaining sections deal strictly with second order theory of station- 
ary random sequences. The spectral representation is of particular interest because 
extensions of this theory are used in the next chapter to define the joint spectral 
density of a pair of random sequences. Many excellent references are available that 
cover this material, so well known theorems are stated without proof except where a 
brief sketch will provide insight to the theory without diverging too far from the flow 
of the discussion [8, 9, 10, 11, 1, 7, 12], 

2.1 Distribution Functions 

Consider a complex valued random sequence {*(*)},* = 0, ±1, ±2, • • • ± ( n - 1) 
Since (X(t)} is complex valued, each element of the sequence can be written as 

= U(t) 4 - j i 

where {U(t)} and {U(f)} are real valued random variables. For finite n, {A'(f)} is 
completely characterized by its 2n-dimens,onal probability distribution function, 

Fu ^) ^n).V(f 1 )„.,V(f n )(Ul, . . . , u n , Vi r J = 

P M*') ^ H'n) < «n, V(t x ) V(t n ) < Vn ). 

( 2 . 1 ) 
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where (x 1 ,...,i n ) e C and x, _ u , + jv, [13]. For notational convenience, the 
distribution function for a complex valued {*(«)} will be written as 

• -»*n). (2.2) 

by which is meant (2.1). When it is not necessary to distinguish whether a set of n 
random variables corresponds to n-length sequence of real valued random variables or 

a j -length sequence of complex valued random variables, the notation (Xi, . . . , X n ) 
will be used. 

A function F is defined as a distribution function if 
1. F is monotonically increasing in each of its arguments, 

F 1® right continuous in each of its arguments, 

3- lim x , — oo F(x iy . . . ,x n ) = 0 (1 < t < n), 

4. F(+oo, . . . , -(-oo) < OO, 

0. F(b u ...,b n ) - E?=1 Px + E,<J Pi] =F •• ■ + (— l) n F(ai, • • • , a n ) > 0, 

for any a, and 6, (i = 1 ,...,n), where p tJ ,„ k is F( Ci,...,c„) with c, = a„ Cj = 

a ’'---' Ck = ° k * = h for all / ^ [9]. The left hand side of the 

condition in item (5) above is equal to the probability that (Xi, . . . ,X n ) falls within 
the parallelepiped 


that is, 


a % < X, < 6 t , 


P ( a i < Xi < b u ...,a n < X n < b n ). 

If F(+oo , . . . , +oo) = 1 then F is a probability distribution function. 

The existence of a stochastic process corresponding to a particular a priori 
of finite dimensional distribution functions, 


F *\ ( z i). ^r,,Xa(*i, xa), • • - , F Xl 


,x n (xi,...,x„), 


i* *s 


is assured by Kolmogorov’s theorem if two conditions are met [14]. That is, a set 
of probability distribution functions are the distribution functions of some stochastic 
process if and only if the symmetry (or permutation ) condition, 

Fx ix ,...,x in (x,, , . . . , x,J = F Xu . „ Xn {x u ...,x n ), 

holds for all permutations of (1,. and the compatibility (or consis- 

tency) condition 


Fx ' *■(*!» •••»*«») = F Xl i B (z 1 ,...,x fn ,+oo,... J +oo), 

holds for m < n. The compatibility condition ensures that the lower dimensional 
distributions coincide with the margins of the n-dimensional distribution. 

2.2 Density Functions 

The derivative of the distribution function, if it exists, gives the density function 


/(*!,•• - ,X n ) 


dFjxu^ • -,x n ) 
dx j . . . dx n 


Explicit reference to the random variables {X,} has been dropped from the notation 

for simplicity when the meaning is clear from the context. A density function has the 
following properties [9]: 


1- /(xj,...,x n ) > 0. 

2. P(ai < Xi < bj , . . . , a„ < X„ < b n ) = / 6 " • 

3 - S-o o' • ■ f^°oof( x u ■ ..,X n )dx l ... dx n = 1. 



• -dx n . 


As a consequence of the compatibility condition 


- J /(xj, . . . , x n )dx m+ i ■ ■ • dx n , 

j —OO J — X> 

SO that the marginal density of (X„ .... X„ } is obtained by integrating out the 
dom variables {X m+I X n }. 


ran 


2.3 Characteristic Functions 

The characteristic function of a random variable X is defined as 

*<") = = J_ x e Juz dF(x), (2 . 3 ) 

where F(x) is the distribution function of X. 

Theorem 2.1 (Bochner-Khinchin) A function <&(w,,..., w „) with $(0 0) = 

1 is a characteristic function if and only if it is non-negative definite. 

In other words, the Bochner-Khinchin theorem says that the class of characteristic 

functions coincides with the class of non-negative definite functions, which by defini- 
tion means that 

n n 

EE $ K-«o^;>o, ( 9 4 ) 

j=i k=i v ~' ; 

for any integer n, where are vectors in W and c„...,c n are arbitrary 

complex numbers. Any function that satisfies (2.4) is a valid characteristic function 
for some random variable. Using (2.4) it is easy to show that non-negative definite 
functions have the following properties: 

1 . $( 0 ,..., 0 ) > 0 , 

2. u, n ), 

3- *(0,...,0)>|#( Wii ..., W w)|. 

A unique inverse relationship exists between a characteristic function $(u^, . 
and its corresponding distribution function F{ i { . . , . , x n ) given by [9], 


p ( a i < 3Ci <b\, ...,a n <X n < b n ) = 


lim 

T — oo 


1 f T fT " 

(27r) n J-T 


_ e -jw fc 6 fc 
J'+'k 


$(^1, . . . ,o>„) <Lo\ . . . (Lfr, . 


2.4 Stieltjes Integral 

The integral on the right-hand side of (2.3) is the so-called Stieltjes integral which 
will be used extensively in what follows [9], The Stieltjes integral of a function g(x) 
with respect to the distribution function F(x) is defined as 

where 


a — x 0 < xj < • • • < x^r = b, and x*_i < x fc < Xk . 

This generalization of the integral remains well defined even when F(x) is not every- 
where differentiable. The improper Stieltjes integral where the interval of integration 
goes from — oo to -)-oo is defined in the usual way as 

lim 

a— •> — o 
6 — ►+<> 

If F(x) is differentiable with respect to x, then the Stieltjes integral reduces to the 
ordinary integral 

/ g(x)f{x)dx 

Ja 

where /(x) = ^1. The Stieltjes integral can be applied to either Riemann or 
Lebesgue integration in which case it is called the Riemann-Stieltjes or Lebesgue- 
Stieltjes integral respectively. 


fg{x)dF(x). 


2.5 Moments 

It is often more convenient to work with a partial characterization of a random 
sequence in terms of a finite set of statistical moments. The first moment is simpU 
the mean and is defined as 

li(t) = E[X(t)}= I" xdF X(t) {x). 

J — OO 
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The second moment is defined as 


U 


R(h,t 2 ) - E[X*{t x )X{t 2 )] = f_ oo f_ oo Xix 2 dF X(timta) ( XuXi ) i (2.5) 

and is called the autocorrelation function. Alternatively, the second central moment 

can be formed by taking the autocorrelation of {*(<) - n(t)} in which case it is called 
the autocovariance, 


C(ti,t2) = R(t u t 2 ) - El^E^ti)}. 

The convention of applying the term autocorrelation to (2.5) is standard among the 
engineering community, however; mathematicians use the term autocorrelation to 
signify the normalized autocovariance. Henceforth, use of the term autocorrelation 
always refers to (2.5), and the mean of {X(t)} is assumed to be zero so that the 
autocorrelation and autocovariance are identical. This assumption results in no loss 
of generality. 


2.6 Stationarity 

A random sequence, (A'(f)}, is strict-sense stationary if its distribution functions 
are independent of shifts in the index t 

Fx ^) ••>*»»») = Rx(t l+ ,).....x(t m +,)(xi, . . . , x m ), 

for any integer i [4). If the joint distribution of two sequences {X(t),Y{t)} is sta- 
tionary then the sequences are said to be jointly (or mutually) stationary. Clearly 
from (2.1) a complex sequence is stationary if its real and imaginary parts are jointly 
stationary. 

A less restrictive form of stationarity, called wide-sense stationarity, only requires 
that the mean be constant 


ONioww, ; 

OF POOR 


1 A > 


{( | • 


M‘) = £[.Y(l)] = /i, 


and the correlation function depend only on the difference of the indexes 

R(U,h)= R(t) = E{X'(t,)X( t, + r)|, 

where r = f, - Consequently, the autocorrelation is independent of the absolute 
starting point, Henceforth, when a random sequence is called stationary it means 
stationary in the wide-sense. Since Gaussian random sequences are fully described in 

terms of their first two moments, in this special case wide-sense stationarity implies 
strict-sense stationarity. 

2.7 Properties of Autocorrelation Functions 

Theorem 2.2 A function R{r) defined on integers is the autocorrelation function of 

a zero-mean, stationary random sequence if and only if it is non-negative definite. 
That is, 

E E R U - Qcjci > o, 

;=l fc=l 

for any positive integer n and arbitrary complex numbers c* [18]. 

As a consequence of being a non-negative definite function, the autocorrelation of a 
complex valued random sequence has the following properties: 

1. R( 0) > 0, 

2- l*(r)| < R(0), 

3. R(r) = RT(-t). 

From the first and second property, the autocorrelation must have its maximum valur 

at Zero la « be non-negative at that point. The third property indicates that Hi - 
is a Hermitian function. 

Comparing the properties of autocorrelation functions with the properties i 
characteristic functions from Section 2.3, it is apparent that an close relation^.,;, 
exists between the two. A variation of Herglotz’ Theorem (Theorem 2.3 in the UN- 
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section) states that a function, R(t), defined on integers is non-negative definite, and 
thus an autocorrelation function, if and only if it coincides on the set of integers with 
a characteristic function, $(w), defined by 

*(«) = [ W e ] “ x dF(x), 

J — TT 

where u> is continuous on R 1 [15], 

2.8 Spectral Representation of Correlation Functions 

Theorem 2.3 (Herglotz) A complex valued function R(t) is non-negative definite 
if and only if it has a representation 

= //'*” dFW ' (2.6) 

where F(\) has the properties of a distribution function on (— 7r,x], 

F(\) is called the spectral distribution function and as stated in the theorem it has 
all the usual properties of a distribution function except that it is defined on (-7 r, jr], 
such that 


F{ A) = 0, A < -7T, 

F( A) = F( tt), A > 7t. 

If a correlation function is absolutely summable, 

£ \R(r)\ 2 <oo, 

T-— OO 

then R(t ) can be represented in terms of it Fourier coefficients 

R(r) = j\^f(\)d\, 

where 

!W = A E c-’' r M(T), 

“ r = -00 

is the power spectral density which has all the properties of a density function defined 

On ( — 7T, tt] . 

i H > > •: . 

CT : if-* -* ' i'*. 1 »2T 


2.9 Spectral Representation of Stationary Random Sequences 

Theorem 2.4 (Cramer) Any zero mean stationary random sequence can be repre- 
sented in the form 

*W = /VVZ ( A), (2.7) 

where (Z( A)} is a complex valued random process with orthogonal increments defined 
on (—7 r, 7rJ . 

A proof of this theorem due to Cramer is given below in a non-rigorous form to 
help illustrate the essential properties of the spectral representation [8, 10, 1, 16, 17] 

Thorough treatments of the Hilbert space theory necessary to understand what follows 
are available in the literature [12, 17]. 

Let H be the Hilbert space formed by the collection of all complex valued random 
variables with zero mean and finite variance. Then for each value of t, the random 
variable X(t) belongs to H. Let H x denote the closed linear subspace of H that is 

spanned by the random sequence {*(<)}. The inner product of two elements of H z 
is defined as 

(*i,*2) = E[x‘ x x 2 \, (2.8) 

and the distance between two elements is defined as 

d(xi t Xi) = ||*i — x 2 || = yj ~E [|xi — x 2 | 2 ]. (2.9) 

Denote by L 2 (F) the set of all complex valued functions, 0(A), on the interval 
( — n i T \ for which the Lebesgue-Stieltjes integral 

y_J0( A)| 2 dF(A), 

exists and is finite. Here F is the spectral distribution function of X(t). Then L 2 ( F i 
forms a Hilbert space Hg with the inner product of two elements 9 X and 0 2 defined 

(0u0 2 )= [* 

J — n 


and the distance between 9 X and 9 2 defined as 

d(0i,0 2 ) = ||*i -0 2 1| = [/j0x(A) -0 2 (A)| 2 dF(A)J’ . 

A linear, one-to-one mapping between two spaces that preserves inner products 
is known as a congruence. A congruence between H x and H e can be established as 

follows. Let *(<) € H t and 6 t { A) = e^‘ € H, be corresponding elements in the two 
spaces. By Herglotz’ Theorem, 

R xx (t,u) = E[x*(t)x(u)] = f e- jXt e~ jXu dF( A). 

Comparing the spectral representation of the autocorrelation with the definition of 

the inner product in H x and H e it is evident that the inner product is preserved for 
a mapping M such that 


M [e"J = i(t), 

and for linear combinations of {0 t (A)} 


M 






Now let 1 (*,,a,](A) be the indicator function defined such that 


Simple functions in 
following manner: 


l(»,.A,lW = 

He can be written 


< A ^ Aj 
otherwise. 

in terms of the indicator 


function 


in the 


~ (2. ID . 

t=l 

where A,_, < A < A, and -t = A 0 < A t < • • • < A n = tt. The indicator function its.-if 
is an element of H$. Define Z( A) as the random process in H x corresponding to th.- 
indicator function in H g through the mapping M as follows: 


Z(\ 2 )-Z(\ 1 ) = A/[1 (Ai>Aj1 (A)]. 
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Clearly, Z( A) is orthogonal over disjoint increments since for [A x , A 2 ) D [A 3 , A 4 ) = 0, 

E[{Z( A 4 ) - Z(X 3 )}'{Z(X 2 ) - Z( A x )}] = <1 ( a 3> a 4 ](A), 1 ( a 1 ,a j] (A)) 

= 0 . 

If inner product is preserved by the mapping M, that implies that distance is also 
preserved so that, 

£ [|Z(A,) - Z(A,)| 2 ) = |IW,,W 

f \ 2 

= / dF(X) 

J ' ' 1 

= F(A 2 )-F(A,). 

Applying the mapping M to (2.10) gives 

A/[0 t (A)] = Jim X]^(A)A/[1 ( a,_ 1 ,a 1 ,(A)] = Jim ^0 ( (A) (Z(A,_ 1 ) - Z(A,)] 

1=0 t=0 

which converges to the Stieltjes integral 

*(o = r o t (\)dz(\). 

J — TT 

Substituting 9 t ( A) = t jM gives the spectral representation in (2.7). The integral on 
the right-hand side of (2.7) is a stochastic Stieltjes integral with respect to a random 
measure, in this case Z(X). 

2.10 Orthogonal Increment Process 
An orthogonal increment process, Z( A), has the following properties: 

1. E[Z( A)] = 0, 

2. Var[Z( A)] < oc. 

3. Cov[Z(X 4 ) — Z(Xj), Z(A 2 ) — Z{\\)\ = 0, (2.11) 

where (A lt A 2 ] D (A 3 , A 4 ] = 0 and C<n;[.Y, Y] and Var[X ] are defined as 

Cov[X,Y\ = E[X m Y ], 

sxmir^i gsp Var[X ] = Cov[X,X]. 
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The orthogonal increment process, Z( A), is related to the spectral distribution func- 
tion in (2.6), F( A), by 

F(\ 2 ) - F( Aj) = Var[Z{ A 2 ) - Z{\i)). (2.12) 

Recall that F(—n) = 0, so letting Aj = — ir causes (2.12) to reduce to 

F(\) = Var[Z{\)-Z{-T)\, 

where Z(—tt) is the starting value of the orthogonal increment process which can also 
be set to zero without affecting the values of the increments [12], Considering (2.11) 
and (2.12), the relationship between Z( A) and F{ A) can be expressed as 

E[dZ-{\)dZ{\)\ = S( A - A)dF(A). 

The inversion formula for Z( A) is given by 

1 JL — (\ V 

Z ( A > = 2r If Zjp-W) + 

It can be shown that X(t ) is Gaussian if and only if Z( A) - Z(-tt) is Gaussian [8], 

2.11 Brownian Motion 

One dimensional Brownian motion [12], {5(A), -*■ < A < ir}, is a simple 
example of an orthogonal increment process which is Gaussian distributed with 

W1 = 0 

Var[B(X)\ = T > . 

The spectral distribution function corresponding to a Brownian motion process spec 
ified as above is given by 


ro 


F(A) = 


2 l±±il 

2ir 


A < —IT, 

— T < A < 7T, 
A > 7T. 
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Substituting F( X) into (2.6) gives the autocorrelation 

R( r ) = r e ]Xr dF{ A), 

J —7C 

a 2 r* .. 

= cr 2 6(h). 


as 


As expected, the autocorrelation of a random sequence with a Brownian orthogonal 

increments process is an impulse at zero lag because the spectral density is uniform 
on ( — 7T, 7rJ. 


CHAPTER 3 

JOINT SPECTRAL DENSITY 


In Chapter 2, a univariate random sequence was shown to have a spectral rep- 
resentation in terms of an orthogonal increments process. Furthermore, taking the 
variance of this orthogonal increments process gives the spectral distribution func- 
tion which acts like a probability distribution function of the frequency content of 
the random sequence. At times it is more convenient to work with the distribution 
function rather than the density function. This should cause no difficulty since the 
density function can be found by talcing the derivative of the distribution function 
(assuming that it is differentiable). 

The purpose of this chapter is to extend the idea of a PDF of frequency content 
to the case of two random sequences. Clearly if the joint spectral density is to be 
completely analogous to a joint PDF, it must be a real valued function of multiple 
arguments (in this case 2). In Section 3.1 classical spectral analysis of bivariate 
random sequences is discussed in terms of the correlation matrix and the spectral 
density matrix. Since these are matrix valued functions of a single argument, they 
cannot be considered candidates for the JSD. And yet, they contain all the available 
information about the joint probability structure of the two sequences. In Section 3.2 
the spectral representation of stationary random fields is discussed and it is noted 
that the spectral density function of a stationary random field has the same form and 
properties that are desired for the JSD. The properties of two-dimensional orthogonal 
increments processes are explored in Section 3.3. The joint spectral density is formally 
introduced in Section 3.4 by treating the autospectra of the bivariate random sequence 
as the marginal spectra of a stationary random field. 
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3.1 Bivariate Correlation and Spectral Properties 
The correlation theory of bivariate random sequences is well documented in the 
literature [1, 12, 7, 19]. Consider two jointly stationary random sequences (A’i(i)},{A’ 2 (<)} 
that have the following spectral representations: 


X, (t) = (' e’ M dZ,(\), 

J — 7T 

(3.1) 

Xi(t) = f e‘ M dZ,( A), 

J —IT 

(3.2) 

where (Zi(A)}, (Z 2 (A)} are orthogonal increment processes. Then assuming that 
both process have zero mean (i.e., E[X x (t)\ = £'[A' 2 (<)] = 0), the autocorrelation 
functions can be written as follows: 

Rn(r) = + r)] = f e> Xr dF n ( A), 

J —IT 

(3.3) 

R 27 (t) = E[X;(t)X 2 (t + r)] = r e^ T dF 22 ( A), 

J — 7T 

(3.4) 


where Fn(A) and F 22 (A) are the spectral distributions of the two processes. Recall 
from Section 2.9 that 


Var[dZi(A)] = dF u { A) 

(3.5) 

Var[dZ 2 { A)] = dF 22 (A). 

(3.6) 

A facet of the analysis of multivariate processes that is absent when considering 
the univariate case is the interaction between pairs of processes. This information is 
captured by the cross-correlation. The cross-correlation is given by 

Ri 2 (t,u) = E[X;(t)X 2 (u)}. 

(3.7 

Substituting (3.1) and (3.2) into (3.7) gives 


Ru(t,u) = f J* e~ }Xt e j ' Xu E[dZ;(\)dZ 2 {\)\. 

(3 ' 
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If {A'i(i)}, {X 2 (t)} are jointly 2” d -order stationary then their cross-correlation will be 
a function only of the shift index r = u-t and will be independent of the starting 
index t. For this to be the case in (3.8), the following condition must hold: 

E[dZ'{\)dZ 2 (X )] = 0 A ^ A, ( 3 . 9 ) 

and therefore, 


Rn(r) = f e’ Xr E[dZ;(\)dZ 2 {\)} = T e^dF l2 { A), (3.10) 

IT J — 


where F\ 2 (A) is the cross-spectral distribution function. The cross-correlation function 
is not necessarily non-negative definite, so it need not share the same properties as 
the autocorrelation. One significant difference is that the cross-correlation may have 
its maximum value away from zero lag. As a result, F 12 (A) is in general complex 
valued and even if it is real valued is not necessarily non-decreasing. It is however a 
function of bounded variation because its real and imaginary parts are the difference 
of two non-decreasing functions respectively [ 7 ], 

For bivariate random sequences, the auto- and cross-correlations are often as- 


sembled into a correlation matrix given by 

. . , -ftn(T-) 

R(r) = 


(3.11) 


Rh{t) R22(t) 

The two cross-correlation elements of the correlation matrix yield redundant infor- 
mation since 


Rl2(r) = 


If the elements of R(r) are absolutely summable. 

£ |£i;(r)| 2 < 00 , i,j = 1,2, 

T = -00 

then spectral density matrix is given by 


f(A) = — y 

{ } 27T T f^ 


-j\t 


R(t) 


(3.12) 


(3.13) 


(3.14) 
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where 


and 


f(A) = 


/11(A) / 12 ( A) 

/21(A) /22(A) 


/«i(A) = 


dF tJ ( A) 
dA 


Properties of the elements of the spectral density matrix include 


(3.15) 


(3.16) 


/11(A) > 0, 

/22(A) > 0, 

I/12I 2 < /n(A)/ 22 (A), 

where the third relation is a consequence of Schwartz’ inequality. The cross-spectral 
density, /12, is also complex valued in general, so it can be written in terms of pairs of 
real valued functions in a couple of different ways. Each of these ways of writing /i 2 
offers a particular interpretive viewpoint. Breaking / i2 down into its real and imag- 
inary parts gives the co-spectrum (real part) and quadrature spectrum (imaginary 
part) as 


/12(A) — 012(A) —7912(A). 

Writing /12 in polar form gives the cross-amplitude spectrum and phase spectrum as 

/i2(A)=a, 2 (A)e J *”< A \ 

where 

“12(A) = |/, 2 (A)| = i/ci 2 (A) + < 7 i 2 (A), 

/ 912(A) \ 

V c«(A) ) • 


<£i2(A) = arg(/i 2 (A)) = arctan 
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The properties and interpretation of the co-spectrum, quadrature spectrum, cross- 
amplitude spectrum, and the phase spectrum are treated in depth in the literature 
[1, 7, 19]. 

The complex coherence is defined as 

7„(A) = -■■*»<*> . 

The complex coherence is sometimes referred to as the correlation coefficient in the 
frequency domain, but it should not be confused with a correlation coefficient between 
the frequency content in the two random sequences. It is, rather, “the correlation 
coefficient between the random coefficients of the components” in the two random 
sequences at each frequency [1]. 

The correlation matrix and the spectral density matrix describe the joint prob- 
ability structure of the bivariate random sequence, however; if the auto-spectra are 
interpreted as the probability densities of the frequency in each random sequence, 
then it might be possible to reorganize the information available in these matrix val- 
ued functions to form a joint probability density function of the frequency content in 
the two sequences. In the following two sections this problem is set within the frame- 
work of the spectral representation of stationary random fields and their marginal 
spectra. 

3.2 Stationary Random Fields 

It is well known that the spectral representation theorem for random sequences 
can also be extended to stationary random fields [20, 7, 1]. A random field is called 
stationary if the mean is constant 

= E[X{t,u)\ = fi, (3.17) 

and the correlation function is function only of shifts in its arguments 


R(t,u,v,w) = R(t,v) = E[X*(t, u)X(v, to)], 


(3.18) 
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where r — v — t and v = w — u. This two-dimensional autocorrelation function has 
a spectral representation 


R(t,v)= r [* e*l Xir+X3V) dF( A x ,A 2 ), 

J — IT J — IT 


(3.19) 


where F(A!,A 2 ) has the properties of a distribution function on (— x, 7 r] x (— ir,ic]. 
If F(Xi,\ 2 ) is differentiable, then the random field has a spectral density, /(Ai,A 2 ), 
given by 


/(Ai,A 2 ) 


dF(A l7 A 2 ) 

dX\dX 2 ' 


where /(A!,A 2 ) has the properties of a density function on (— x, it] x (— 7 r, 7 r]. 


Theorem 3.1 Any zero mean stationary random field can be represented in the form 

X(t, u) = r r ^ x ' t+x ^dZ( A x ,A 2 ), (3.20) 

J —7 T J — 7T 

where {Z( A l7 A 2 )} is an orthogonal increment process on ( — 7r, 7 t] x (— 7 r, 7 r] [20, 7j. 


The proof of this theorem given below is the author’s extension of Cramer’s proof of 
Theorem 2.7 to random fields. 

Let H x be the closed linear subspace of H that is spanned by the random field 
{^f(t,u)} (see Section 2.9). The inner product and distance for H x are defined as 
(2.8) and (2.9) respectively. 

Denote by L 2 (G) the set of all complex valued functions, ^(Ai,A 2 ), on ( — 7r,7r] x 
( — 7r, 7r] for which the two-dimensional Lebesgue-Stieltjes integral 

Jjy^X^dGi^X^ (3.21) 

exists and is finite. Here G is the two-dimensional spectral distribution function of 
X(t,u). Then L 2 (G) forms a Hilbert space H w with the inner product of elements 
and x[ 2 defined as 

(01,0a) = I* r 0,*(A 1 ,A 2/ V 2 (A 1 ,A 2 )<fG(A 1 ,A 2 ). 

J _7r J — 7f 


(3.22) 



To construct a congruence between H x and H let X(t , u) € H z and ^ (U (A 1 , A 2 ) = 
e j(*i <+a 2 «) g be corresponding elements in the two Hilbert spaces. By Herglotz 
theorem in two dimensions, 

Rx*(tyU,V,W) = 

E[x'(t,u)x(v,w)\ = r r e-il Xit+x > u) e ]{XlV+x > w) dG(\ u \ 2 ) 

J — IT J —IT 

— f f 0* u (Ai , A 2 )V’vu/(Ax, AjJd^Aj, A 2 ) 

J — 7T 7 — IT 

{x(t,u),x’(v,w)) = (l/>tu,Tpvw)- 

Then the inner product is preserved for the mapping T such that 

T J e i(A,t+A 2 «)j _ X (t,u). 

Let A 2 ) be the indicator function defined on rectangles such that 


) A 2 ) — 


1 p, < A| < fij, Vi < A 2 < Vj 
0 otherwise. 


Simple functions can be constructed from this two-dimensional version of the 
indicator function as follows: 


lMAi,A 2 ) = Jhn J2 V’»*(Ai,A2)l( w . 1 , M .]x(^_ If ^](Ai,A 2 ) 

«,i=i 

where 


fti-i <Ai < A 4 *, 
i'i-l <A 2 < i/;, 


r Aq ^ ... ^ A n — tt , 

— 7T = Vq < ■ • • < V n = 7T. 


The indicator function is also an element of H \p, so we can define its correspondine 
random process, Z(Aj,A 2 ), in H x as 


- Z(n 2 ,v 1 ) - Z(fi U v 2 ) + Z(nx, v x ) = T [l( /llf/1J jx(i/ ll ^](Ai, A 2 )] .(3. J ! 



It can be shown that Z( Ai, A 2 ) has orthogonal increments since for non-intersecting 
rectangles, the inner product of the indicator function is equal to zero. Since distance 
is preserved by the mapping T, 


E 


| Z(fi 2 ,i/ 2 ) - Z{n 2 ,v 1 ) - Z(ni,i/ 2 ) + Z(/z 1 ,t/ 1 )| 2 

^2)|| 

r^i /m 2 

= / / dG{ \ u \ 2 ) 

= G(n 2 , v 2 ) — G(fi 2 ) ^ 1 ) — G(ni, v 2 ) + G(fx 1 , i^i). 
Applying the mapping T to (3.23) gives 
T A 2 )] = Jiin ^2 A 2 )T 1(^,_, iMi ]x(v 7 _,,i/ ; ](Ai, A 2 )J , 

<,j=x 

n 

= JiiS, 21 M^u a 2 ) [Z(m, vj) - Z{m, u^) - z{^ u Vj ) + 

>,i=i 

which converges to the two-dimensional Stieltjes integral 

x(t,u)= r r MXuX'mxuw- 

J — 7T •/ — IT 

Letting 0tu(Ai, A 2 ) = gives the spectral representation of a two-dimensional 

random field 


X(t,u)= f r e>l x ' t+x " ) dZ( A x ,A 2 ). 

*/ —7T J — 7T 


(3.24) 


3.3 Properties of Z(Ai,A 2 ) 

An orthogonal increment process in two dimensions has properties similar to the 
one dimensional case. 


1 . £;[z(a 1 ,a 2 )] = o, 

2. Var[Z(A 1 , A 2 )] < 00 , 

3. Cov [AZ{(/i 3 ,^ 4 ] x (t/ 3 ,i/ 4 ]} ,AZ {(/*!, /z 2 ] x (i/x,t/ 2 ]}] = 0, 
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where (^i,/i 2 ] x (iy u i/ 2 ] D (/x 3 ,^ 4 ] x (v 3,1/4] = 0 . The increment operator A in two 
dimensions is defined as follows: 

AZ {(/il , /i 2 ] X (j / l,i / 2 ]} Z(fJ, 2 , vf) Z(flx, l/ 2 ) — Z(fl 2 , V i ) ~h Z(/Ji, l/i). 

The random measure AZ{(pi, p 2 \ x (i/i,i/ 2 ]} is related to the spectral measure by 

A F{(mi,M 2} x (v x,v 2 ]) = Var[AZ{{p u p 2 ] x (*/,,**]}]. 

To see this relationship more clearly, consider computing 

Cov [AZ {(p 3 ,p 4 J x [v 3 ,v 4 ]} , A Z {(p 1 ,p 2 J x (i/i,i/ 2 ]}] , 

where (^1 , ^ 2 ] x , i/ 2 ] H (^3,^4] x (1/3, 1/4] ^ 0 . This situation is depicted in Fig- 
ure 3 . 1 . The only non-zero contribution to the covariance is from the region where 
the rectangular increments overlap (the shaded region in the figure). Therefore, the 
covariance reduces to 

Cov [AZ {(/i 3 ,/i 4 ] x (1/3, i/ 4 ]} , A Z {(/ii,/i 2 ] x (i/i,i/ 2 ]}] = Var[AZ {{p 3 ,p 2 \ x (v\, v A }}\ . 

As the rectangular increments are shrunk down to infinitesimals, the rectangles will 
only overlap if their vertices coincide so that we are left with the short-hand notation. 

Cov[dZ{ Aj, A 2 ), dZ( A*i, A 2 )] = < 5 (Ai — Aj)« 5 (A 2 — \ 2 )dF(\i, A 2 ). 

3.4 Joint Spectral Density 

Definition 3.1 The joint spectral distribution, F{ Ai,A 2 ), of two random sequences. 
Xx(t) and X 2 (t) with corresponding orthogonal increment process Z\(X) and Z 2 { A) is 
equal to the spectral distribution function of some stationary two-dimensional random 
field such that 

F( X,7r) = Var[Z l (X)), 

F(ir,X) = Var[Z 2 (\)], 

and is consistent with the Cov[Zi(X\), Z 2 (X 2 )] oway from the margins. 
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ement process. 




31 


The last line of the definition deserves some comment. In general, the joint spectral 
distribution (or density) is underdetermined by X x (t) and X 2 (t). However, under 
certain conditions it is possible to establish a relationship between the cross spectral 
distribution, Fi 2 (A), and F(X x ,X 2 ) at points away from the margins. In these cases, 
the joint spectral distribution can be uniquely determined for a given pair of random 
sequences. Initially it might seem that this idea violates the well known maxim 
that the joint probability density can not be uniquely determined from the marginal 
densities. However, there is more information available than what is contained in the 
marginal spectra. Recall that the power spectral density is not a sufficient statistic for 
a random sequence since the phase information is discarded by taking the magnitude 
squared of the orthogonal increment process, Z( A). Thus conditions may exist under 
which the cross-correlation or equivalently the cross-spectrum contains the additional 
information necessary to form the joint density. 

To construct the joint spectral density of a bivariate random sequence, it is 
necessary to take the information in the spectral density matrix, (3.15), and somehow 
form it into a two-dimensional spectral density function like /(Ax,A 2 ). Consider the 
spectral representation of the autocorrelation of some random field, X(t,u), that was 
given in (3.19) 

r(t,v)= r r^^dF(x u x 2) . 

J — r J — ir 

Notice that setting either r or v equal to zero is equivalent to integrating out Ax or 
A 2 respectively, 

Rn(r) = R(r,0) = fJ\ }X ' T dF( A,,A 2 ) = J\> x ' r dF u ( A x ), (3.25) 

R 22 (t) = R(0,t) = T f e jX * T dF(Xi, A 2 ) = r e jXiT dF 22 (X 2 ), (3.26) 

J —7 T J — 7T J — IT 

where f?n(r) and R 22 (t) are the marginal correlation functions and f r ’n(A 1 ) and 
F 22 ( A 2 ) are the marginal spectral distributions. This is necessary for T’(Ai,A 2 ) to 
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Figure 3.2. Example of a two-dimensional distribution function, F(Ai, A 2 ). 
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satisfy the compatibility condition. Then Fu(\\) and F 2 2 ( A 2 ) are related to F(Ai, A 2 ) 

by 

F ll (X) = Var[Z l (X)-Z 1 (-ir)) 

= Var[Z( A, 7r) — Z( A, —n) — Z(—7r,7r) + Z( — 7r,—7r)] 

— F(X, 7r), (3.27) 

F 22 (X) = Var[Z 2 (X) - Z 2 (-n)) 

= Var[Z(7r, A) — Z(w, — tt) — Z(— 7r, A) + Z(— tt, — 7r)] 

= F(tt,A). (3.28) 

If the marginal spectral distributions of a stationary random field are set equal to 
the auto-spectral distributions of a pair of random sequences, then F( Ai,A 2 ) will be 
consistent with the joint spectral distribution of those sequences along the boundaries 
Aj = 7r and A 2 = -jt. Figure 3.2 shows a plot of a two-dimensional distribution function 
on (— 7r,7r] x (— 7r,7r]. Notice that on the boundaries along Ai = — tt and A 2 = — tt 
the joint distribution function is equal to zero and that it is monotonically increasing 
in each of its arguments until they reach 7r where the joint distribution function is 
equal to the marginal distribution function. To complete the specification of the joint 
spectral distribution, it remains to determine F(Aj,A 2 ) in the interior region away 
from the margins in terms of the joint statistics of X\{t) and X 2 (t). 

Under certain symmetry conditions, knowing the diagonal slice, F( A, A), in addi- 
tion to the margins is sufficient to completely specify F(Aj, A 2 ) over ( — 7r,7r] x ( — 7r, ttJ. 
The Gaussian distribution is an example since its density function exhibits elliptical 
symmetry (see Appendix B). The projection-slice theorem is an interesting property 
that relates slices through a two-dimensional correlation function to projections in 
its Fourier transform as depicted in Figure 3.3 [21]. If R(r,v) is evaluated along the 
slice formed by setting v = r, then 

R(t) = R(t,t) = r r e* x ' +x ' )r dF(\ u \ 2 ), (3.29) 

J — tt J — n 
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results. Taking the one-dimensional Fourier transform of R(r ) yields a projection of 


the joint power spectral density along the 45° line formed by setting A 2 = A^ 

1 00 - 

/(A) = 5= £ mr)'-» r (3.30) 

Z7r -00 

= ^E e ~ jXT [/" // j(A,+Aj)t ^(Ai,A 2 )] . (3.31) 

Assuming .F(Ai,A 2 ) is differentiable, i.e., 

dF(\\, A 2 ) = /(Ai , A 2 ) d\i d\ 2 , (3.32) 

then 

1 00 r r* /*» 

/(A) = —Y, e _,ir [£ £ « J,A,+X!lr /(Ai, A 2 ) dA, dA 2 

= s££[|^ ,+ ^ , 1^„A l)<i A,dA 2 

= ^ (3.33) 

If /(At,A 2 ) is separable in its arguments such that 

/(A i,A 2 ) = / 1 (A 1 )/ 2 (A 2 ), 


then (3.33) reduces to the convolution of £(Ai) and / 2 (A 2 ). Notice that in this case. 
/(A) corresponds to the probability density function of the sum of Ai and A 2 where 
they are statistically independent random variables. 

The author knows of no existing method for computing arbitrary slices R(t,v) 
that are not on the margins. At this point we reach an impasse and need to find 
a way to further constrain the problem. The interrelationships among the various 
time and spectral domain representations are summarized in Figure 3.4. Lines that 
have arrows on both ends signify Fourier type operations which are bidirectional. 
Lines with arrows on just one end signify one-way operations that are either of the 
covariance type, or involve extracting the margins. Starting from X\ (t) and X 2 (t) 
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at the bottom notice that there is no direct path to the goal of a joint spectral 
distribution at the top. The premise is that given certain restrictions on X\ (t) and 
Xi (<), the combined information in the auto- and cross-spectra (or correlations) can 
form a link to F( Aj^) (or R(t, u)). 







3 


i » A, 2 ) Z(X.i,X2) 



Figure 3.4. Schematic diagram of the relationships among various time-domain and 
spectral representations. 



CHAPTER 4 


MARGINAL SEQUENCES 

4.1 Introduction 

One way to constrain a pair of random sequences so that they uniquely specify 
the joint spectra density is to assume that they are the marginal sequences of some 
stationary random field. This restriction implies that the spectra of the two sequences 
correspond to the marginal spectra of the random field, but in general the converse is 
not true. That is, if the spectra of two sequences are the marginal spectra of a station- 
ary random field, it does not necessarily follow that the sequences are the marginal 
sequences of that or any other stationary random field. This is a consequence of the 
“many-to-one” relationship between sequences and their power spectral densities. It 
is well known that many sequences may share the same power spectral density, but 
that any particular sequence has one and only one power spectral density [22]. 

4.2 Basic Concepts 

The marginal sequences of a two-dimensional stationary random field, {X{t, u)}, 
are defined as the pair of random sequences that result from setting t = 0 and u = 0 
respectively, 


X x {t) = X{t, 0) 

X 2 (u) = .Y(0,u). 

Figure 4.1 shows how the marginal process are related to the random field as a whole. 
One of the distinguishing properties of marginal processes that is clearly illustrated 
in Figure 4.1 is that the two sequences must start at the same value. 


X 1 (0) = A' 2 (0) = A(0,0). 
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Figure 4.1. Marginal process of a two-dimensional random field. 


Further restrictions result directly from the assumption that the random field is 
stationary. 

Theorem 4.1 The marginal sequences of a stationary random field are individually 
stationary. 

This theorem follows directly from the definitions of a stationary random field and 
its marginal processes. A random field, {A^(f,u)} is stationary iff 

f?n(r,u) = E[X*(t, u)X(t + r,u + u)]. (4.1 ) 

Letting u = 0 and v = 0, gives 

R X x(t, 0) = E[X m (t,Q)X(t + t, 0)], 


= J F[A7(0A 1 (< + r)], 

= Rn(r). 


(4.2) 
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Similarly it can be shown that X 2 (t) is also stationary. From (4.2) it is clear that the 
autocorrelations of the marginal processes correspond to the margins of the autocor- 
relation of the stationary random field. That is, 

Ru(t) = Rxx(t, 0) 

R2i{t) = Rxx(0,t). 

Theorem 4.2 The marginal sequences of a stationary random field, X(t,u), are 
jointly stationary iff the autocorrelation of X(t,u) can be written in the form 

Rxx{t,v) = Rxxiv - r). 

Again, the proof follows directly from the definitions of a stationary random field and 
its marginal processes. In (4.1) let u = 0 and t + r = 0, then 

i2xx(r,v) = J E[X*(-T,0)AT(0,u)], 

= £[*,’(-r)Jf a (t;)], 

= R, 2 (-T,v). (4.3) 

Only in the special case where R xx {r, v) = Rx X {v — r) will the marginal sequences be 
jointly stationary. Returning to (4.1), if we let t = 0 and u + v = 0, then 

R xx (t,v) = E{X’( 0,-u)X(r,0)], 

= E[Xl(-v)X, (r)], 

= R 21 (-v,t). (4.4) 

Comparing (4.3) and (4.4), gives 

Rn{~T,v) = R 2 i(-v,t). 

Since R xx (t, v) = R’ xx {—r , — u), we also have 

Ri2(~t,v) = Kiffr, -u) 

R2\(~v,t) = R^i{v, -r). 
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Therefore R 12 and R 2 1 are hermitian functions and are non-negative definite. As a 
result, F x2 and F 2X are real-valued and non- decreasing. 

4.3 Properties in the Spectral Domain 

It might be expected that the correlation relationships for marginal processes 
from the previous section have a significant effect on the spectral domain characteri- 
zations. In this section that is shown to indeed be that case. 

Theorem 4.3 If Xi(t), i = 1,2 are the marginal sequences of a stationary random 
field, X(t,u), then 


Z\( A) = Z( A, 7r) - Z( A, -tt), (4.5) 

Z 2 {\) — Z(tt,A) — Z(— tt, A), (4.6) 

where Z,( A) are the orthogonal increment processes corresponding to the marginal 
sequences and Z(\ x ,\ 2 ) is the orthogonal increment process for X(t,u). 

To prove this theorem, consider the spectral representation of a stationary ran- 
dom field, 

X(t,u) = /l /l ej(A,, + A?U> ^( Al ^2)- 

Setting u = 0 gives an expression for the marginal sequence in terms of the two- 
dimensional orthogonal increment process, Z(Aj.A 2 ), 

X,(t) = X(t, 0) = „A 2 ), 

= <<Z ( A1 ’ A *), < 4 - 7 ) 

where the inner integral is with respect to \ 2 and outer integral is with respect to Aj. 
Recalling that the spectral representation of -Vi(f) is given by 

*i(<) = f/ XldZ x{ A), 
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and comparing this with (4.7) gives a relationship between the marginal orthogo- 
nal increment process, Zj(A) and the two-dimensional orthogonal increment process, 
A 2 ), 


dZi( Xi)= r dZ{ A lt A 2 ), 

J—ir 


where the integral on the right hand side is with respect to A 2 . This expression can 
be evaluated by integrating both sides with respect to Aj over the interval (a, 6] and 
applying the definition Stieltjes integral. 

[ b dz l (\ l )= f b r dz( a 1 ; a 2 ) 

Ja J a J — ir 

n 

Z\{b) - Zi(a) = lim Z (^ v i) ~ Z(ia,vj- 1 ) - Z{m- U uj) + Vj - X ) 

where a = jx 0 < • ■ • < fi n = b and — tt = i/ 0 < • • • < v n = tt. By the noting the 
cancelation of terms in the summation, it is easy to verify that 


Zi(b) — Z-i(a) = Z(b, ir) — Z(b, —i r) — Z(a , i r) + Z(a, — t). 

Equating terms on both sides gives 

Zi(b) = Z(6,7r) — Z(b, — 7r) 

Zi(a) = Z(a,Tr) — Z(a, — ir). 

In a similar manner (4.6) can be proved. 

Corollary 4.4 If Z\{\) and Z 2 ( A) are orthogonal increment processes corresponding 
to the random sequences X\ (t) and X 2 (f) respectively, and are related to the two 
dimensional orthogonal increment process, Z(Aj,A 2 ), by 


Z 1 (A) = Z(A,7r)-Z(A,-7r), 
Z 2 (A) = Z(7r,A)-Z(-rr,A), 



then 


X l (t) = X(t, 0), 
X 2 (t) = X(0,t), 
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where X{t,u) is the stationary random field corresponding to Z(Ai,A 2 ). 

The proof is a simple variation on the proof of Theorem 4.3 and will not be given 
here. 

Setting the starting point of the orthogonal increment processes to zero which 
results in no loss of generality, we get the following expressions: 

Z,( A) = Z(A.t) 

Z 2 ( A) = Z(;r,A). 

Substituting Z( A, 7 r) and Z(tt, A) into the expression for the cross-spectral distribution 
gives 


f 12 (A) = E[z;{\)z 2 {\)\ 

= Cou[Z(A,7r),Z(7r, A)]. 

Recall that Z(Aj, A 2 ) is a two-dimensional orthogonal increment process so the covari- 
ance of Z(A l5 A 2 ) over two intervals is only non-zero where the two intervals overlap. 
If the two intervals are disjoint then the covariance is identically zero. Therefore, 
referring to Figure 4.2, it is easy to show that 

Cov[Z(\ U Tr), Z(rr, A 2 )] = F(\ u A 2 ). (4.8) 


Hence for this special case, the cross spectral distribution is equal to the joint spectral 
distribution along the diagonal Aj = A 2 . 
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Figure 4.2. Domain of a two-dimensional orthogonal increment process. 


4.4 Jointly Stationary Marginal Sequences 
To evaluate the effect that assuming the marginal sequences are jointly stationary 
has on the JSD, start with the expression for the correlation of a stationary random 


field 


R(t, v) = E u)X(t -f- r, u + u)] . 

Suppose R(t,v) is evaluated along v = 0, then 

/?(r,0) = £[r(f ) uP + r,u)], 

= £[X*(0,0)X(r,0)], 

= £[ x ;( o )*,( r )] = e[x;(o)Xi(t)). 


fln(r) 


«2l(r) 


(4.9) 
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Similarly, if R(t,v) is evaluated along r = 0, then 

R( 0,v) = £(X*(t,u)X(t,u + v)], 

= E[X'( 0,0)Jf(0,«)J, 

= £( v ,-( 0 ) a : 2 ( u )1 = E{x;m)X 2 (v)]. ( 4 . 10 ) 

V ' ' 

^12(v) H22 (u) 

Also since both Xi(0) and .A 2 (0) are equal to AT(0,0), 

Rn(r) = E [A7(0)AT 1 (t)] = ^(O^r)] = R 22 (t). (4.11) 

Comparing (4.11) to (4.9) and (4.10) it is recognized that all of the auto- and cross- 
correlations functions are identical, 


^n( r ) — Ru( t ) = R2i(t) = /?22( r )- 


In the spectral domain this implies that 


Fi i(A) = F n ( A) = F 21 ( A) = F 22 ( A). 


Recalling that 


F 12 (A) = F(A,A), 


we have 


F(A,7r) = F(A,A) = F(7r,A). 

Unfortunately this is a trivial result from the standpoint of estimating a joint spectral 
density given an arbitrary pair of random sequences because the two sequences are 
required to have the same PSD. Even so, there is some pedagogical benefit to carrying 
this a bit further with a couple of simple examples. 
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4.5 Some Examples 

A class of two-dimensional distribution functions defined on (— 7r,7r] x ( — 7r, tt] 
satisfying the conditions, 


F( A,tt) = F( A, A) = Ff7r,A), 


can be constructed as follows: 


1. Choose any valid one-dimensional distribution function (7(A), 

2. Set F(A 1? A 2 ) = 6 r (min(A 1 , A 2 )). 

Since F(Ai, A 2 ) must be non-decreasing in each of its arguments, it will have contours 
of constant value as shown in Figure 4.3. The joint spectral density is given by 


/( Ai,A 2 ) 


F(Ai,A 2 ) 
d\\ dX 2 

G(min(At, A 2 )) 
dX\ dX 2 


Note that the min(A x ,A 2 ) can be written as 


min(Ai,A 2 ) — -(Ai + ^2 l^i — ^ 2 !)- 


By applying the theory of generalized functions (also known as the theory of distri- 
butions) we can evaluate the derivative of min(Aj,A 2 ) with respect to Ai and A 2 as 
follows: 


-£^E-min(Ai,A 2 ) = + A, - |A, - A,|), 

= + sgntAt - A 2 )), 

= r-TT-( 1 + (2u(A x — A 2 ) — 1)), 

L d\\ 

= -1 U( A, - A 2 >, 


= 6(X t -A 2 ), 


where sgn(x) is the signum function and u( x) is the unit step function [23, 24, 25]. 
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Figure 4.3. Contour plot of F , (A 1 ,A 2 ) for jointly stationary marginal sequences. 


Applying the chain rule of differentiation, it can be shown that 
/(A i, A 2 ) = ^(min(A 1 , A 2 ))£(A! - A 2 ), 

where ^(A) is the density function corresponding to (7(A). The joint spectral density 
is concentrated on the line Ai = A 2 and is zero everywhere else. The following two 
examples illustrate this class of JSD. 

4.5.1 Uniform Margins 

The distribution function, .F(A 1 ,A 2 ), when the marginal spectra are uniform i' 
given by 


F( Ai,A 2 ) 


< 


0 

1 (At + A 2 — \\i — A 2 |) + ttJ 

<7 2 (A]+t) 

2ir 

g 2 ( A; + t) 

2 IT 


Ai or A 2 < — tt, 

—ir < Ai, A 2 < 7 r, 

— IT < At < 7T, A 2 > 7T, 
Ai > 7 T, — 7T < A < 7 r, 
\i and A 2 > 7 r. 
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This function is plotted in Figure 4.4. The marginal spectral distributions F u (^i) 
and ^ 22 (^ 2 ) 3X6 equal to each other and correspond to a standard uniform distribution 
on (— 7r,7r]. The joint spectral distribution, however, is not strictly uniform in the 
sense that the joint spectral density is not flat over ( — 7r, 7 t] x ( — 7r, tt] . 

In order to find the autocorrelation function R(t.v), first evaluate 

dFj A X ,A 2 ) 
dX 1 d A 2 
a 2 6 ( Aj - A 2 ) 

2tt 

The spectral density corresponds to a uniform density concentrated along the diagonal 
Ai = A 2 . 


/(A,,A 2 ) 


Substituting /( Ai,A 2 ) into (3.19) gives 

R(r,v) = ^ ~ X *) dX ' dX >' 

= ^~ r e j(T+v)A2 d\ 2 . 

27T J — ir 


Therefore, 



e j(r+v)A 2 ^ 


2ir t — —v 

< 

0 otherwise. 


R(t, v) = o- 2 6(t 4 - v). 

The autocorrelation is thus a train of impulses running along the diagonal r = — v. 


4.5.2 Gaussian Margins 

The distribution function, F(A!,A 2 ), for the case of Gaussian marginal spectra 
with zero mean and unit variance is given by 


Q 


+ H ^ I 

<7 


Q(^) 

a 


Q(^) 


1 


Ai or A 2 < — 7 r, 

— 7T < A lf A 2 < 7T, 


— 7T < Ai < 7T, A 2 > 7 T, 
Aj > 7T, — 7T < A < 7T, 
A! and A 2 > tt. 


F(X u X 2 ) = { 
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Figure 4.4. Plot of F( Ai,A 2 ) for jointly stationary marginal sequences with uniform 
marginal spectra. 
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See Appendix B for the definition of the (^-function. 

The marginal spectral distributions F n (A 1 ) and F 2 2 (A 2 ) are equal to each other 
and correspond to a univariate Gaussian distribution. Figure 4.5 shows a plot of this 
spectral distribution function. 

The joint spectral density is given by 

dF(X u X 2 ) 


/(A,,A 2 ) = 


d\idX 2 

1 (*1 -»)* 

e <5(Ai — A 2 ). 


\/2t' 


tct 


Substituting into (3.19) gives 

R(t,v) = -tL- r r e j l x ' T+x ' v) e- L ^ 1 8(X 1 - X 2 )dX,dX 2 , 

\/ Z7TCT J-* 

= -jL- r 4x„ 

V 2 ircr J-t 

-S e -5* 2 ( T + v ) 2 c Mr+t^ 

The envelope of this autocorrelation function is shown in Figure 4.6. 

4.6 Jointly Harmonizable Marginal Sequences 
Assuming that the marginal sequences are jointly stationary places an unaccept- 
able constraint on the joint spectral density. This constraint can be relaxed however 
by considering a generalization of stationarity called harmonizability [15]. A harmo- 
nizable random sequence can be represented by 

X(t) = r e> Xt dZ( A), 

J — n 

but whereas in the stationary case, Z( A) has orthogonal increments, for harmonizable 
sequences, the increments are not necessarily orthogonal. The autocorrelation is given 
by 

R(t,u) = E[X m (t)X(u)] 

e- ]Xt e jXu E[dZ'(X)dZ{X )] 


fir rn 

J — ir J — ir 
= r r e j{Xu ~ Xt) dF{ A, A), 

J — w j — 7 T 


(4.12) 



lambdal 


1 


-1 


lambda? 




3^-3 


Figure 4.5. Plot of F{ Ai,A 2 ) for jointly stationary marginal sequences with Gaussian 
marginal spectra. 
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Figure 4.6. Plot of R(t, v) for jointly stationary marginal sequences with Gaussian 
marginal spectra. 
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where .F(A,A) is of bounded variation and in general complex valued. In the special 
case where X(t) is stationary, (4.12) reduces to (2.6). 

Suppose X\(t) and X^it) are the marginal sequences of a stationary random 
field X(t,u) and they are individually stationary, but jointly harmonizable. Then 
the autocorrelations will be one parameter functions of shifts in t, but the cross- 
correlations will be two parameter functions. Consider once again, the autocorrelation 
of the stationary random field 

R(t, v) = E [X“(£, u)X{t 4- r, u + u)] . 

Letting u = 0 and t = — r, gives 

R(t,v) = E [X'( — t, 0)Af(0, u)] , 

= E[X{{-t)X 2 { 0,v)], 

= Ru(-t,v). 

Similarly, letting t = 0 and u = —v, gives 

R(t,v) = E[X‘(0,v)X(t,0)], 

= E[X;(-v)X l (r)}, 

= R2l(-V, t). 

Note that the cross- correlations are non-stationary, so R\i(t, u) and /? 2 i(r, u) can 
NOT be reduced to Ru{v — r) and Rn(v — r) respectively. 

Recalling that the autocorrelation of a stationary random field is a Hermitian 
function i.e., 


R(t,v) = R’( — t, -u), 


it is easy to show that 


R\2{~T,V) - R- U (T. -V) 
R2i(-v,t) = /^(u, -r). 
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Therefore, the cross-spectral distributions are guaranteed to be real valued and non- 
decreasing. 

4.7 Some Examples 


4.7.1 Generalized Brownian Motion 

Perhaps the simplest example of a two-dimensional orthogonal increment process 
is a variation on the Brownian sheet on ( — 7r,7 r] x (—7 r, x]. Let {Z(X\, A 2 ), -t r < 
Ai,A 2 < 7T } be an orthogonal increment process that is Gaussian distributed, with 


E[Z(X U X 2 )] = 0, 


Var[Z( Aj, A 2 )] — — - (Ai + 7 t)(A 2 + r) + j-j (A 2 — i r 2 ) (A 2 — 7r 2 ) 


47T 2 


(4.13) 


(4.14) 


Then the corresponding distribution function, F(Aj,A 2 ), is given by 


^(Ai,A 2 )= { 


0 Aj or A 2 < — 7r, 

[( Ai -I- 7r)(A 2 + 7r ) + 4 ^T (A 2 — 7T 2 ) (A 2 — 7T 2 )J — T < A X , A 2 < 7T , 

— 7T < Aj < 7T, A 2 > 7T, 
Aj > 7T, — 7T < A < 7 IT, 
cr* Aj and A 2 > 7r. 


2?r 

^(Aa+y) 

2?r 

r2 


The marginal spectral distributions Fi^Ax) and -^ 22 (^ 2 ) are equal to each other 
and correspond to a uniform marginal spectral density. In general, the joint spectral 
distribution is not strictly uniform in the sense that the joint spectral density is not flat 
over (— 7 r, 7 r] x ( — 7r,7r] unless the parameter k is set to zero. This particular family 
of generalized bivariate uniform density functions is called Morgenstern’s bivariate 
uniform distribution (see Appendix A). It should be noted that other families of 
bivariate density functions exist that generate uniform marginal densities [26]. 

The joint spectral density is given by 

dF( \ u \ 2 ) 


/( Ai,A 2 ) = 


dX\ dA 2 

r 2 





00 


Substituting into (3.19) gives 
2 

r ( t ^) = £^ F_Jl/ (XiT+x ^ ] ( 1 + ^ t AiA2 ) dXidX2 ' 

= £ £ £ eJ,A,T+A "’ "*"* + 0 £ £ ".ix. 

2 

= <7 2 ^(r,u) + ^p; J_ *ie jXlT d\ij X 2 e jX * v dX 2 . (4.15) 


The two integrals on the right hand side of (4.15) evaluate to 

f Xe^ T dX = - 4 r [sin(r7r) — T7rcos(r7r)] 
J-TT T l 


(4.16) 


A plot of this function for continuous r is shown in Figure 4.7 If r is restricted to 
integers, (4.16) reduces to 


P Xe jXT dX = -^cos(rTr) 

J —ir T 

27r£ 


-(-1) T+1 - 


(4.17) 


A plot of this function is shown in Figure 4.8. 
Substituting (4.16) back into (4.15) gives 


R(t, v) = ct 2 6(t, v ) — 


cr 2 n 


47r 4 r 2 i; 2 

<7 2 K 


sin(r7r) — T7rcos(r7r)] [sin(t>7r) — u7rcos(u7r)] 


= c 2 &{t, v ) [sin(r7r) sin(t»7r) — tit cos(t7t) sin(u7r) — 

4 7T 4 T j V^ *■ 


vtt sin(r7r) cos(ti7r) + tvtt 2 cos(ttt) cos(l>7t)| . 


Restricting r and v to integers, reduces this to 


R(t,v) = <7 2 6{t,v) - (-1) 


T + V 


ct 2 k 

47 t 2 tv 


(4.18) 


Evaluating slices in R(t , v) corresponding to t = 0 and v = 0 gives 


Rn(r) = R{t, 0 ) = & 2 6(t), 
R 22 (v) = R(0,v) = cr 2 6(v), 








and the cross-correlation is given by 


■) i 


Ru(r,v) = R(-r,v) = <t 2 S(t,v ) -f (-1) 


T+V 


<J K 


Ak 2 tv 

The two-dimensional autocorrelation for the generalized Brownian sheet is shown in 
Figure 4.9 with the exception of the impulse function at the origin. The parameter 
k is proportional to the correlation coefficient between Aj and A 2 . Notice that as 
k — ► 0, the oscillations away from zero lag vanish and the autocorrelation becomes 
an impulse at the origin. This is to be expected because for k = 0 the JSD is flat 
over ( — 7T,7r] X ( — 7T, 7r]. 


4.7.2 Bivariate Gaussian 

In the generalized Brownian motion example, even though the marginal se- 
quences are not jointly stationary, the marginal spectra are identical. Let’s consider 
a case where the marginal spectra are different from each other. If the joint spectral 
density is bivariate Gaussian then the marginal spectra will be univariate Gaussian, 
but in general they will have different means and variances. In the following example 
the means are set to zero for simplicity, but the variances are unequal. 

Let {Z( Ai,A 2 ),— 7r < Aj,A 2 < 7r } be an orthogonal increment process that is 
Gaussian distributed, with 


E[Z{ A 1 ,A 2 )] = 0, 

yar[Z(A lt A 2 )]:=z(^Ap) +q(~) + Q (-) - L (-, —,p) - 1. 

V<Tl <7 2 J \(TiJ \cr 2 ) \<7! cr 2 J 

The corresponding spectral distribution, F(Aj,A 2 ), is given by 

r 

0 Ai or A 2 < —7 r, 

L (£, %,p) + Q (£) +Q(*)-L £,/>) - 1 ^ A„ A, < T. 

F(A„ A,) = { L (k i,,) + Q (i) + Q (i) - L (i, £,„) - 1 —ir < A, < x, A . 

L fe. £./>)+« fe) - Q (£) - L (£• *r) ~ 1 

1 Ai and A 2 > 7r. 
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Figure 4.9. Plot of R(t, v) for continuous r and v. 


59 


See Appendix B for the definitions of the Q-function and the T-function. Unlike the 
uniform distribution, the bivariate Gaussian distribution is non-zero over all of R 2 , 
so the tails of the Gaussian will be assumed to be sufficiently small at the boundaries 
of (— 7T, 7rJ x (— 7T, tt] that 


«(- 

Wi 


<7i/ 

<7 2 ) 


) *0, 


Q (-) w o, 

\<7 2 J 


% * \ 

,01 cr 2 ) 


LI^Ap 

\cr i cr 2 


0, 

% 0, 

TT TT \ 

— , — ,pj « 0 . 

<?\ a 2 ) 

Otherwise the truncated bivariate Gaussian distribution should be used [27]. If this 
assumption is valid, then i r (Ai,A 2 ) reduces to 

0 Ai or A 2 < — 7r, 

t t5’p) + Q (?fj + Q (^) — i — 


F( A.,A 2 ) = | 




—TT < A i < 7T, A 2 > 7T, 
Al > 7 r, — TT < A < TT, 
Aj and A 2 > tt. 


The joint spectral density is given by 


A? 2pX\X 2 X\ 


L 2 ( 1 _ / > 2 )V cr i a i a 2 


An 

°l) 


/(Al ’ Ai) = 2^,Wl-^ CXP 

Substituting into (3.19) gives 

R( t » y ) = exp [— i(cxfT 2 + 2pa x cr 2 Tv + <rfu 2 )j 
The autocorrelations of the marginal sequences are given by 


— TT < Ai, A 2 < T 


Rn(r) = R(t, 0) = e-W T \ 
R22 (t) = R(0,t ) = e-H^, 
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and the cross-correlation is given by 

Ru{t, v) = R{—t, v) = exp | — ^(<7j r 2 - 2/7cr,cr 2 ru -f- a\v 2 ) 

See Appendix B for plots of the bivariate Gaussian. 

4.8 Simulating Marginal Spectra with Particular JSDs 
Marginal sequences with specified joint spectral densities can be simulated using 
a variation on a technique proposed by Zrnic for generating I&Q sequences with 
weatherlike spectra [28]. Writing the complex valued random field in terms of the 
inphase and quadrature (I&Q) components, 

X(t, u) = I(t,u) +jQ(t,u), 

they can be modeled by 

I(t,u) = s(t,u) cos <f>(t,u) + n(t, u) cos xp(t, u) (4.19) 

Q(t,u) = s(t, u) sin <b(t, u) + n(t, u) sin u ) (4.20) 

where s(t,u) and n(t,u) axe Rayleigh distributed signal and noise envelopes and 
4>{t,u) and 0(t,u) are uniformly distributed phases on [0, 2x], The noise envelope, 
n(£,it), is assumed to be broadband compared to s(t,u). Expressing the I&Q fields 
in terms of a two dimensional Discrete Fourier Series (DFS) gives 

i M m 

/(«,») +><?(*,«) = T75 E E \/ P ( k . ') (4.21) 

where P{k,l) is the instantaneous power of the signal plus noise and 8(k,l) is a 
uniformly distributed phase on [0, 2tt\. If the shape of the true JSD is given by S(k,l) 
and the noise power per discrete frequency is a constant, N, then P(k,l ) is given by 

P(k, l) = - [S(k, l) + N] In U(k, /), (4.22) 

where U(k, l) is a uniform random number on [0, 1]. The steps for generating the l&Q 
field can be summarized as follows: 
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1. Select a desired spectral shape S(k,l). 

2. Set the signal-to-noise ratio (SNR). 

3. Take the inverse discrete Fourier transform of the product of \JP(k, l) and an 
independent phase e^ k ' l K 

The marginal sequences are easily extracted from the resulting I&Q field by setting 

Xi{t) = X(t, 0), 

x 2 (t) = x( 0 , 0 - 

In the following two subsections, this simulation technique is used to generate the 
marginal sequences for the examples in subsections 4.5.2 and 4.7.2 respectively. 

4.8.1 Jointly Stationary with Gaussian Marginal Spectra 
In this example, the marginal sequences are jointly stationary with Gaussian 
marginal spectra. The ideal shape of the joint spectral density is 

1 (*i-«i 2 

S( A 1? A 2 ) = -—e—te-6{ X l - A 2 ). 

yzTca 

Figure 4.10 shows the dB magnitude of the true JSD of the signal plus noise for a SNR 
of lOOdB. The quadratic shape of the log of the Gaussian can be clearly distinguished 
rising out of the noise floor along \\ = A 2 . 

Since the marginal sequences are jointly stationary, the auto-spectra are identical 
and equal to 

fn( A) = / 22 (A) = - X 2 ). 

V 2xa 

A plot of the auto-spectra of the simulated marginal sequences is shown in Figure 4. 1 1 . 
The joint spectral density is formed by setting /(A X ,A 2 ) equal to /n(A) along the 
diagonal and zero elsewhere. 

If the assumption that the marginal sequences are jointly stationary is dropped, 
the non-stationary cross-correlation can be computed by taking the ensemble average 
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lambda 


Figure 4 . 11 . Plot of the marginal spectra /11(A) and /22(A) 







(4.23) 


over N independent realizations of the marginal sequences (29] 

*.<!..« )*«<*.)• 

t 

Figure 4.12 shows the estimate of the JSD based on the two dimensional Fourier 
transform of the non-stationary cross-correlation where the ensemble average is taken 
over N = 1000 realizations. A ridge along Aj = A 2 is still visible, but it is embedded 
in a structure dictated by the shape of the marginal spectra. For a single realization 
(N = 1), the non-stationary cross spectral estimate is given by 

/l2(Ai,A 2 ) = /ll(^l)/22(A 2 ). 

The structure characteristic of these estimates of the non-stationary cross-correlation 
is an artifact of this outer-product of the sample auto-spectra. In the limit as N — ► oc, 
/( Ai , A 2 ) should approach the true JSD, but the convergence does not appear to be 
particularly fast. 

4.8.2 Jointly Harmonizable with Bivariate Gaussian JSD 
The case where the marginal sequences are jointly stationary is not a realistic one 
for estimating the JSD from the non-stationary cross-correlation since the JSD can be 
found much easier from the auto-spectra. Consider a more practical example where 
the marginal sequences are jointly harmonizable with a bivariate Gaussian shaped 
JSD. For this example the ideal shape of the joint spectral density is 

1 / A? 2pAiA 2 A| \ 

2(1 — P 2 ) “ <7K7 2 + A J] ' 

In the simulations, crj = <r 2 and p = .95, so that the presence of a correlation in the 
joint spectrum would be clearly visible. Figures 4.13 and 4.14 show a 3-D mesh plot 
and a contour plot respectively of the dB magnitude of the true JSD of the signal 
plus noise for a SNR of lOOdB. 

Using (4.23) to estimate the JSD. the plot shown in Figure 4.15 was generated. 
It is more difficult to see, but the desired gaussian shape is embedded in the central 
mound. It can be seen more clearly in the contour plot of Figure 4.16. 


5(A 1 ,A 2 ) = 


27T<Ti<T 2 v / 1 — p 2 


exp 
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Figure 4.12. Plot of the JSD, /( Ai,A 2 ). 
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Figure 4.15. Plot of /(Ai,A 2 ). 
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4.9 Summary 

In this chapter, marginal sequences are shown to be a class of sequences for which 
the JSD can be determined. Two drawbacks are immediately apparent. First, for 
the two sequences to be valid marginals, they must start at the same value. It might 
be tempting to try to relax this condition by setting the starting values of the two 
sequences to zero, or the average of the two starting values, or employing some other 
replacement scheme, but that would violate the condition that the two sequences be 
individual stationarity. Further research needs to be done applying these sequences 
to practical applications. The other drawback is that jointly stationary marginal 
sequences necessarily have identical PSDs. The resulting JSD is a degenerate case. 
A more general class of marginal sequences can be treated by assuming them to be 
jointly harmonizable. Then estimating the JSD involves estimating the non-stationary 
cross-correlation which is quite difficult in practice. In the simulations, it was assumed 
that an arbitrary number of realizations of the marginal sequences were available for 
computing ensemble averages. 



CHAPTER 5 


AN EXPONENTIAL MODEL 

An alternative approach to constraining the sequences is to specify a signal model 
that explicitly incorporates the joint spectral density in its definition. Consider a pair 
of random sequences {X(t),Y(t)} that can be modeled by 

*(0 = £<*.e j(A ‘‘ + * ,) + 6, (5.1) 

1=1 

= +Cu (5.2) 

1=1 

where {a,} and {/?,} are constant complex valued amplitudes, {A,} and {i/,} are 
random frequencies, {#,} and {<^,} are independent uniform random phases, and f t 
and Ct are zero mean independent white Gaussian noise sequences. 

Using this model, a link can be established between the bivariate random se- 
quence and a family of random fields for which the spectral representation is equal 
to the JSD of {X(t), Y(t)}. For this purpose, it is necessary to assume that pairs of 
frequencies {A,, r/,} are independent and identically distributed according to a bivari- 
ate probability density function, p Xu ( A,i/), that coincides with the normalized joint 
spectral density, f Xv { A,i/). Let the characteristic function of p Xu (\,v) be denoted by 
Marginal functions for /a„(A,i/), p Xu { A.t/), and v) are written in the 

usual manner. 

A similar model is often encountered in the familiar parameter estimation prob- 
lem in which devising “good” estimators for the amplitudes, frequencies, or phases 
is of interest [30]. This model is different in one significant respect. The frequency 
content in (5.1) and (5.2) is a random variable, and the goal here is not to estimate 
the frequency (or amplitude, or phase). Rather, it is shown below that if a pair of 



random sequences can be modeled by (5.1) and (5.2), then it is possible, under certain 
conditions, to determine the joint probability density of A and v. 


5.1 Correlation and Spectral Properties 
The first question to answer is whether X(t) and Y(t) are individually and/or 
jointly stationary. The mean of X(t) is given by 

K 


E[X(t)\ = E 


2> 




+ it 


Li=l 


= X>£ [e j(A * t+<? ')] + £[6], 

t = l 

= [e jA ' £ ] E [e*] 


i=i 
= 0. 


The last step results from the fact that E = 0 for 0, uniformly distributed on 

(— 7 r, 7 r]. Similarly, it is easy to show that E\Y{t)\ = 0. The autocorrelation of Af(<) 
is given by 


R xx (t,u) = E[X'(t)X(u )], 


= E 


K 


X>: 




L \t=i 
K K 


■'«•■) (s 


OLkt' 




+ iu 


= £ £ a i ( a ^ E [e" jA, ‘e jAtU ] E + N^(u - t), 


i=l i=fc 

where is the mean white noise power for £< . Since the phases are independent and 
uniformly distributed, 


E 


e : (*k-0.) 


1 k = i, 

0 otherwise, 


and the autocorrelation reduces to 
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where r = t — s and therefore X(t) is wide-sense stationary. Notice that the expec- 
tation on the right-hand side of (5.3) is, by definition, the characteristic function of 
A, [4]. Since the {A;} are independent and identically distributed, the characteristic 
function can be pulled outside the summation and the autocorrelation becomes 

R xx (t) = A$\(t ) + N ( 6(t), 


and similarly, 


R vy (r) = + N c 6(t), 


where A = £,=i |q,| 2 , B = £,=i |/?,| 2 , and jV c is the white noise power of Ct- 
The cross-correlation of X(t) and Y(t) is given by 


R xy (t,u) = E[X-(t)Y(u)], 


(jt ct’e- i ^ t+8 ' ) + z'j /W ( *' fcU+ * fc) + C u j 


= E 

K K 

= J2Y< a iPkE 

i=i *=1 


e -jA,t e 


£ e A<h c-*i) 


(5.4) 


If 0, and <f> k are uncorrelated for ail z, k then the cross-correlation is zero. If, on the 
other hand, 0 t - and 4> k are allowed to be correlated for i = k, then the cross-correlation 
will be scaled by a factor determined by the joint probability density of 0, and d> x . For 
example, if the phases axe jointly distributed according to Morgenstern’s bivariate 
uniform (see Appendix A), then 




k - z, 

0 otherwise. 


For our purposes, the constant scale factor can by assumed to be unity without 
loss of generality. This assumption corresponds to setting 0, = <A,. 

Substituting back into (5.4) gives 


K 


R xy {t,u) = 


=J (-A,t+i/,u) 


i=l 
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Note that in general, the sequences are not jointly stationary for A, ^ v t . The expec- 
tation on the right-hand side can be recognized as the joint characteristic function 
of A, and Since the pairs {A,, i/,} are independent and identically distributed, the 
characteristic function can be pulled outside the summation to give 

R xy (t,u) = Gt’ 

If $Ai/(— t,u) = $Ai/(u — t) then the sequences are jointly stationary. 

It is remarkable that the cross-correlation may or may not contain information 
about the joint probability structure of A and v depending on a fairly subtle detail of 
the signal model. Thus it is difficult to make general statements about a relationship 
between the correlation between sequences in the time domain and the correlation 
in their frequency content. As in the case of marginal sequences, if X(t ) and Y(t) 
are jointly stationary, then the joint spectral density is concentrated along A = i/. In 
order to treat a more general class of problems the non-stationary cross-correlation 
must be used. 

Taking the discrete time Fourier transform of Rx X {t), i2 yv (r), and R xy (t, u ) gives 
the auto-spectra as 

= Ahm + 

/„( A) = BfM + 

and the cross-spectral density as 

/xy(A, v) = /a«/( — A, v)- 

Considering that the sequences are composed of exponentials with additive white 
noise, a question arises of whether spectra are discrete or continuous. Since the fre- 
quency content is a random variable, there is no reason to expect that the A,’s (or v x s) 
will be harmonically related and thus the sequences will not be periodic. Nevertheless, 
they are “almost” periodic and taken at face value, a single realization of X{t) (or 



K(f)) would be expected to have a discrete spectrum [1]. However, in the expected 
value, the true spectrum is continuous and equal to the probability density function 
of the frequency. This is a model for which the sample spectrum can be significantly 
different from the true spectrum. For small K, the sample spectrum consists of a 
flat noise component with peaks located at the frequencies of the exponentials. Most 
likely, this will bear little resemblance to the true spectrum. Therefore, in order for 
the sample auto-spectra to look anything like the true spectrum, K must be large. 


5.2 Simulations 

In the following simulations, the joint spectral density is bivariate Gaussian with 
zero mean and equal variances so that 


MKv) = u^i -^ exp 

and the marginal spectra are given by 

fxW = 


2<t 2 (1 -p 2 ) 


— (A 2 - 2p\u + V 2 ) 


27T 


-e 


7TCT 


1 - 2 

fM = 


\/2w< 


7TCT 


Figure 5.1 shows the power spectral density for K = 256 exponentials, and an SNR 
of 50dB. A roughly Gaussian shaped can be distinguished centered at A = |. 

In the next two subsections, estimates of the joint spectral density, /ai/(A, t/), are 
based on estimates of the non-stationary cross-correlation computed as, 


1 


N 


l 

where N is the number of realizations in the ensemble average. Taking the two- 
dimensional FFT of R ry (-t,u) gives /a„(A,i/). The resulting estimates are displayed 


as contour plots so that correlations between A and u will be most visible. 
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Figure 5.1. Power spectral density of X(t) for K=256. 


5.2.1 Single Exponential ( K = 1) 

The simplest case is when X(t) and Y{t) are the sum of a single exponential and 
a white random noise sequence and the amplitudes are set to unity. 

X(t) = e i(A,+s) +6, 

Y{t) = e j{ut+<t,) + Ce- 
The the cross-correlation is given by 

U 6 = <f> then the phases exactly cancel, leaving 

R x y (M) = 4>v^( -a.O- 

Figures 5.2 thru 5.4 show contour plots of /. w ( A. v) for p = 0, 0.5, and 0.9 respectively. 
Ensemble averages are over 1,000 realizations in each case. What should be noticed 
in these plots is that as p increases, the plot become more prolate along A = u. 
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5.2.2 Multiple Exponentials (K >1) 

With more than one exponential in the sum, the situation becomes a bit more 
complicated. Consider the case of two exponentials. 

X(t) = 4- e j{X2t+82) 4- ft, 

Y{t) = e ] ( u ' t+,t> ' ) + e :{u2t+,t,2) + Ct- 

The cross-correlation is given by 

R xy (s, t) = E + e -> (X23+)h) 4- f,) 4- e J ^ t+ ^ + &)] , 

= E 4- E E [e j( ® 2 * 5,) ] + 

E U(^*- Aj,) l E [e j(e,_<?2) ] + E [e j( ^‘- Alj) ] other terms. 

The “other terms” include all the product terms with £, and f t - In the expected value, 
all the terms of the summation go to zero except E e^‘ /,<-A,4 *| and E e J ( l ' 2t-A2 * ) . 
but when estimating the cross-correlation from a finite number of realizations, these 
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“vanishing” terms can contribute a significant error. The more exponentials that are 
included in the model (the higher K is) the more of these cross-product terms appear 
in the estimate. The upshot is that as A' is increased, more realizations are needed 
to obtain a reasonable estimate of R xy (s,t). Figure 5.5 shows the contour plot of 
for K = 2 and p = 0.9 averaged over 1,000 realizations. Notice that in 
Figure 5.5 the contours axe not as prolate as they were for the single exponential case 
shown in Figure 5.4. Many other situations could be considered here, but this is not 
intended to be an exhaustive study of the estimation issues associated with the joint 
spectral density. The insights to be gained from these simulations are that 

1. For this model, it is possible to get the joint spectral density from the non- 
stationary cross-correlation. 

2. A large number of realizations are needed for meaningful estimates. 

5.3 Summary 

The exponential model proposed in this chapter, provides an alternative example 
of how a bivariate random sequence can be constrained so that its joint spectral 
density can be determined. If the phases, 9, and <£,, are uncorrelated then the cross- 
correlation is zero, and can reveal nothing about the joint probability structure of A 
and v. If the phases are correlated, then the cross-correlation is related to the joint 
characteristic function of A and v. Except for some special cases, the cross-correlation 
is non- station ary. When A, = i/ t , the joint spectral density is concentrated along the 
diagonal just as was the case for jointly stationary marginal sequences in Chapter 4. 
In general, simulations showed that the shape of a joint spectral density estimate 
based on the cross-cross correlation, averaged over a sufficient number of realizations, 
does reflect the degree of correlation between A and v. 




Figure 5.5. Contour plot of /ai/(A,^) ( N = 1000, if = 2, p = .9). 



CHAPTER 6 


CONCLUSIONS 
6.1 Motivation 

The original motivation for pursuing the definition of a joint spectral density was 
rooted in a problem from the radar signal processing of weather. In this application, 
random sequences represent the backscatter returns from range cells that extend 
radially from the radar antenna. The frequency content of each random sequence 
is related to windspeed in its corresponding range cell through the Doppler shift 
principle. For detecting hazardous weather conditions with a pulsed Doppler radar, 
the windspeed gradient is often a quantity of particular interest. When estimating a 
windspeed gradient based hazard index, the ability to also estimate the probability 
density of that index would be of considerable value. It is in this context that the 
need for the joint spectral density arose. If, through the joint spectral density, a 
joint probability density of velocity in adjacent range cells can be estimated, then 
random variable transformations can be used to get a PDF for windspeed gradient or 
elementary functions of the windspeed gradient. This method could be implemented 
either parametrically or non-parametrically. If prior knowledge of the form of the 
joint PDF of velocity is available, then its parameters could be estimated from the 
joint spectral density. Otherwise, the joint spectral density could be used directly. 
The contribution of this dissertation has been to develop the theoretical framework 
necessary to define and apply the joint spectral density. 

6.2 Properties of the Joint Spectral Density 
In Chapter 3, the joint spectral density is defined in terms of the two-dimensional 
spectral density of a stationary random field that has marginal spectra corresponding 
to the power spectral densities of the two sequences. There are two stages of u many-to- 
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one’’ relationships operating here. At the first stage, many different random sequences 
can have the same power spectral density. At the second stage, once a two-dimensional 
spectral density has been specified, there are many random fields that can have that 
same power spectral density function. Therefore it is possible to think of the joint 
spectral density as forming a bridge between families of random sequences and families 
of random fields. The spectral theory of multivariate random sequences and the 
spectral theory of random fields axe both well established. The importance of the 
joint spectral density is that it pulls these two bodies of theory together in such a 
way that it is possible to talk about the joint probability density function of the 
frequency content of multiple sequences. Therefore all of the familiar techniques for 
interpreting and manipulating probability density functions can be applied to the 
analysis of the frequency content of multiple random sequences. 

The projection-slice theorem provides a relationship between slices in the two- 
dimensional autocorrelation of this random field, R(t,v), and projections in the joint 
spectral density. For spectral density functions with elliptical symmetry such as the 
bivariate Gaussian, knowing i?(r, 0), R(0, r), and i?(r, r) is sufficient to determine the 
joint spectral density everywhere. The margins of R(t,v) are recognized to be the 
autocorrelations of the bivaxiate random sequence, but further research needs to be 
done to determine if there are any conditions under which R(r,r) can be computed 
from those sequences. 

In the absence of any constraints on the random sequences, the joint spectral 
density is underdetermined. In other words, a gap exits between the conventional 
functions of bivariate spectral analysis, e.g. the auto- and cross-spectral densities, 
and the joint spectral density. In Chapters 4 and 5, two types of sequences are 
identified for which the JSD be determined by realizations of a bivaxiate random 


sequence. 



6.3 Marginal Sequences 

In Chapter 4 it is shown that if the bivariate random sequences are the marginal 
sequences of a stationary random field, then the JSD can be written in terms of the 
cross-spectral density. When the sequences are jointly stationary, it is shown that 
the autospectra of the two sequences are necessarily equal to each other and the JSD 
is concentrated along the diagonal and has the same shape as the autospectra. This 
is considered to be a degenerate case. By relaxing the condition of joint stationarity 
to require only that the sequences be jointly harmonizable, a more general class of 
sequences can be treated. 

Examples are provided of the JSDs for the jointly stationary and jointly harmo- 
nizable cases. A simulation method is presented for generating marginal sequences 
that correspond to an arbitrary joint spectral density specified by the user. For 
a jointly stationary and a jointly harmonizable case, realizations of the simulated 
marginal sequences were used to estimate the JSD. Plots of the true and estimated 
JSD are included for comparison. It is clear from these plots that the estimated 
JSD does not closely resemble the true JSD. However, the effects of a correlation be- 
tween the two frequency variables is visible so the potential for making parameterized 
estimates of the JSD is promising. 

6.4 Exponential Model 

In Chapter 5 the sequences are assumed to be modeled by a sum of complex 
exponentials and additive white Gaussian noise. For a particular sequence, the fre- 
quency of each exponential in the sum is an independent realization of a random 
variable. The random variables that generate the frequencies in the two sequences 
are distributed according to a joint probability density function that corresponds to 
the joint spectral density. If the phases are independent and uniformly distributed 
on (0, 2x], the sequences are shown to be individually stationary and with autocor- 
relations proportional to the marginal characteristic function of the frequency. The 
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cross-correlation is found to be zero if the phases in one sequence are uncorrelated 
with the phases in the other sequence. If the phases are correlated, then the cross- 
correlation is proportional to a reflection of the joint characteristic function of the 
frequencies. 

Simulations are presented for the joint spectral density distributed as a bivariate 
Gaussian. For the case of a single exponential, simulations were run for p = 0, 0.5, 
and 0.9 and contour plots of the estimated joint spectral density reflect the degree 
of correlation between the frequency content of the two sequences. As the number 
of exponentials is increased, it was found that more realizations need to be ensemble 
averaged to mitigate the effects of cross-product terms in the cross-correlation. These 
cross-product terms go to zero in the expected value, but when the expected value 
is estimated by averaging over a finite number of realizations they can introduce a 
significant error. 

6.5 Future Work 

For the marginal sequences and the exponential model, the joint spectral den- 
sity is equal to the non-stationary cross-spectral density reflected about one of the 
frequency variables. As such, estimation of the cross-correlation is a primary con- 
cern. The fact that the cross-correlation is in general non-stationary complicates the 
estimation problem considerably. Since ergodicity does not hold, time averages can 
not be substituted for ensemble averages. This means that, except in certain special 
cases, multiple realizations are needed to get meaningful estimates. A potential rem- 
edy to this estimation problem is the concept of local stationaxity which offers the 
possibility of estimating the non-stationary cross-correlation from a single realization 
[29, 31]. The cross-correlation is called locally stationary if it can be factored into a 
stationary part and a non-stationary part by introducing the substitution 

= t - t 2 = t + -, 



into Rxy(t\ y £ 2 )- Then the cross-correlation can be written as 

Rxy(t u t 2 ) = CT , (^y^) 

where C is a constant, Ti(f) is the instantaneous power, and ^(r) is a stationary 
sequence. If Ti(t) is slowly varying with respect to r 2 (r) then they can be estimated 
separately [29]. This and other implementation issues should be a fruitful area for 
future research. 

6.6 Final Assessments 

In this dissertation a theoretical framework has been developed for considering 
the joint spectral density of bivariate random sequences. A significant contribution 
has been made by defining and exploring the properties of the joint spectral density 
in terms of a link between the spectral representation of bivariate random sequences 
and the spectral representation of stationary random fields. Since the joint spectral 
density was found, in general, to be underdetermined with respect to the two se- 
quences, there is a need to devise constraints on the sequences that allow the JSD 
to be determined. Two such constraints are investigated here, but further research is 
called for to identify other classes of signals that yield unique JSDs. 
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Appendix A 

Families of Bivariate Uniform Distributions 

There is more than one solution to the problem of determining a bivariate den- 
sity function that yields uniform marginal densities. Two one-parameter families of 
bivariate distributions that satisfy the consistency requirement for the marginals are 
the Morgenstern’s uniform distribution and Plackett’s uniform distribution [26]. 

A.l Morgenstern’s Uniform Distribution 
Morgenstern proposed a family of distribution functions given by 

H(x t y) = F(*)G(y){l + k[ 1 - F(x)][l - G(y)]}, -1 < « < 1. (A.l ) 

The corresponding density function is given by 

h(x,y) = f{x)g(y){l + k[2F(x) - l][2G(y) - 1]}, -1 < « < 1. (A. 2) 

If F(x) and G(y) are uniform distributions on [0,1], then Morgenstern’s uniform 
distribution is given by 

h{x,y) = [1 + k(2x- l)(2y- 1)], -1<«<1. (A.3) 


Figures A-l thru A-5 show h(x,y) for various values of k. 


A. 2 Plackett’s Uniform Distribution 

The density function for Plackett’s bivariate uniform distribution is given by 

0[(0 ~ l)(x + V ~ 2 xy) + 1] 


h(x,y) = 


{[1 + (V> - l)(x + y)] 2 - - l)xy}? 


, xl> > 0 . 


(A.l 


Figures A-6 thru A-10 show h(x,y ) for various values of xj>. 
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Figure A-9. Plackett’s uniform density for rf> — 2. 




Figure A- 10. Plackett’s uniform density for = 5. 


Appendix B 

Some Properties of the Gaussian Distribution 
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B.l Univariate Gaussian 
A standard univariate Gaussian is given by 

G(x)=-$=e = A (B.l) 

v27r 

and the distribution function is given by 

F x (x) = P(X < x) = -L f e=£ dt, (B.2) 

V ^ j —GO 

where the integral on the right hand side can not be evaluated in closed form. Instead 
it is often tabulated in terms of the area under the tail of the Gaussian, the so called 
Q function [32] 

1 ro o .2 

Q(x) = -^=jf e-TTdt, (B.3) 

or in terms of the error function 

erf(x) = -^= / e^~ dt. (B.4) 

s/tt Jo 

The relationships among F(x),Q(x), and erf(x) for a standard zero mean, unit vari- 
ance Gaussian function axe as follows: 

F(x) + Q(x) =* 1, F(x) = Q(-x), erf(x) = 2F(v/2x) - 1. 

For a Gaussian random variable y with mean /i and variance er 2 , the distribution 
function, F(y), can be found from the Q and erf functions replacing x with This 
gives 

™=«(^)=H i+ ' rf fe))- (B - 5) 
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B.2 Bivariate Gaussian 
A standard bivariate Gaussian is given by 


g(x,y,p) = 


l 


27rvT^7 J 
and can be written in terms of Z as 


exp 


- (x 1 - 2pxy + V 2 ) 


2(1 — P 2 ) 


(B-6) 


= 7 fh Z(l)Z (7^7) ■ (B ' 7) 

Figures B-l thru B-3 show plots of the bivariate Gaussian for p = 0, p < 0, and 
p > 0 respectively. As in the univariate case, the integral of the bivariate Gaussian is 
tabulated in terms of the area under the tail called the L function [32] 

L(h, k, p) = g(x,y,p)dxdy - Z{x)dx Z(w)dw (B.8) 


where 



Some useful relationships among F, L , and Q are 


F(h, k; p) = F(k, h ; p) = L(-h , -k, p) = L(-k , -h, p ) (B.9) 

= L(h,k,p) + Q(h)+Q(k)-l (B.10) 

F(—h, k,p) = 1 - Q(h) - L(h, k, -p) (B-ll) 

F(h, -k, p) = \ — Q(k) - L(h , k, -p). (B.12) 


B.3 Symmetry Properties of the Bivariate Gaussian 
If a bivariate time series has a Gaussian shaped joint power spectral density, then 
the characteristic function is of the form 

$(n, r 2 ) = e -^ a2T ^ 2pabTxT7+b2T ^ (B.13) 

which has elliptical symmetry resulting from the 

a 2 x 2 + 2 rabxy + b 2 y 2 


( B . 1 4 ) 




Figure B-l. 2-dimensional Gaussian with p = 0 






Figure B-2. 2-dimensional Gaussian with p > 0 



Figure B-3. 2-dimensional Gaussian with p < 0 




form of the argument of the exponential. If r = 0 then (B.14) reduces to 


x 2 v 2 
— + — = c 2 


b 2 


c v 


(B.15) 


which describes a family of ellipses for which the major and minor axes coincide with 
the coordinate axes (see Figure B-4). A non-zero value for r corresponds to a rotation 
of axes by an angle of a, where a is given by 

a 2 -b 2 


COt Q = 


2rab ' 


(B. 16) 


as shown in Figure B-5. 



Figure B-4. Family of ellipses for r = 0. 
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